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ABSTRACT 



Neutron stars in binary orbit emit gravitational waves and spiral slowly together. 
During this inspiral, they are expected to have very little vorticity. It is in fact a 
good approximation to treat the system as having zero vorticity, i.e., as irrotational. 
Because the orbital period is much shorter than the radiation reaction time scale, it is 
also an excellent approximation to treat the system as evolving through a sequence 
of equilibrium states, in each of which the gravitational radiation is neglected. In 
Newtonian gravity, one can simplify the hydrodynamic equations considerably for an 
equilibrium irrotational binary by introducing a velocity potential. The equations 
reduce to a Poisson-like equation for the potential, and a Bernoulli-type integral 
for the density. We show that a similar simplification can be carried out in general 
relativity. The resulting equations are much easier to solve than other formulations of 
the problem. 



Introduction 



Ever since the discovery of the first neutron star binary by Hulse & Taylor (1975), we have 
known that these systems exist and will undergo orbital decay via gravitational wave emission. 
The inspiral and coalescence of binary neutron stars is one of the primary targets for gravitational 



wave detectors now under construction, like LIGO, VIRGO, and GEO (see, e.g., Abramovici et 



al. 1992, or Thorne 1994). And the coalescence of binary neutron stars is the basis for several 



models for gamma-ray bursters (Paczyhski 1986; Eichler et al. 198S; Narayan, Paczyhski, & Piran 



1992). Clearly, a theoretical understanding of coalescing binary neutron stars is an important 



astrophysical problem. 

Analyzing coalescing neutron stars requires the full machinery of general relativity. While 
some issues can be addressed with Newtonian or post-Newtonian gravity, such as the gravitational 
wave form at large separation, others cannot even be posed in these weak-field limits. For 



example, consider the recent controversial claim by Wilson, Mathews, and Marronetti ( Wilson fe 



Mathews 1995| ; |Wilson, Mathews, fc Marronetti 19961 ; |Mathews, Marronetti, fc Wilson 199^) that 
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massive neutron stars collapse to black holes before they merge. Gravitational collapse of a star in 
equilibrium is a consequence of the nonlinear nature of Einstein's equations; it cannot be treated 
correctly with Newtonian gravity. If neutron star binaries do in fact collapse before they merge, 
there are enormous implications for gravitational wave detection: The entire detection strategy is 
based on having accurate templates for the wave forms, and current templates predict relatively 
high sensitivity because without collapse there are many orbital periods before coalescence. And 
of course if the neutron stars collapse, then their final coalescence cannot be a source of gamma 
rays. 

There are other aspects of the inspiral problem that require a fully relativistic treatment. 
For example, if the stars do not collapse before merger, their combined mass is likely to be 
greater than the maximum mass of a cold, rotating neutron star. Then the merged remnant must 
ultimately collapse to a black hole. But it is not clear whether this collapse occurs immediately on 
a dynamical time scale, or whether thermal pressure from shock heating will be sufficient to hold 
the star up for a while. In this case, collapse will occur on a neutrino-dissipation time scale. We 
also do not know how much angular momentum the final remnant will have. What will the fate 
of the system be if the angular momentum exceeds the maximum allowed value for a Kerr black 
hole? Will the excess angular momentum be radiated away gravitationally, or will it be ejected 
in a circumstellar ring of matter, possibly leading to planet formation? All these questions have 
observational implications that can be addressed only by fully relativistic simulations. 

A number of groups have undertaken the construction of general relativistic codes that can 
treat the binary problem. One needs a code that solves Einstein's equations in three spatial 
dimensions plus time, and simultaneously solves the equations of relativistic hydrodynamics. This 
is a formidable challenge, and not surprisingly various approximations to the full problem have 
been proposed. As we have argued above, it seems important not to give up the strong-field 
aspects of the problem by using Newtonian or post-Newtonian gravity. Instead, we will focus 
on two approximations that are likely to be well-satisfied for real neutron star binaries, at least 
up until the orbit becomes unstable and the stars finally plunge together: quasiequilibrium and 
irrotational flow. 

In Newtonian gravity, one can find exact equilibrium states for a binary neutron star system. 
In general relativity, a binary loses energy by gravitational wave emission and the stars spiral 
ever closer together. However, the time scale on which the orbit changes is much longer than the 
orbital period, at least until the orbit becomes unstable at the innermost stable circular orbit. We 
can therefore approximate the evolution as proceeding through a succession of equilibrium states 
of decreasing separation, in each of which gravitational radiation can be neglected. The rate of 
progression along the sequence is determined by the rate of gravitational wave emission, but the 
structure at each point along the sequence can be calculated to excellent approximation ignoring 
the radiation. 

Relativistic quasiequilibrium neutron star binaries have been constructed by Baumgarte et al. 
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(1997a, b,c). For Einstein's equations, they used the approximation scheme of Wilson and Mathews 
([Wilson fc Mathews 1989) ; Wilson 1990| ; [Wilson fc Mathews 1995| ), which is a consistent way 



of neglecting small gravitational radiation terms in Einstein's equations. This scheme has been 
calibrated by Cook, Shapiro, & Teukolsky (1996) on rapidly rotating single neutron stars, where it 
gives excellent results. Wilson and Mathews evolve the binary system by solving the full dynamical 
equations for the matter in the instantaneous background metric, and then updating the metric 
quantities at each time step by re-solving the approximate Einstein's equations. Baumgarte et al., 
on the other hand, take advantage of the quasiequilibrium approximation for the matter equations 
too: They reduce the matter equations to a Bernoulli integral, which can be solved much more 
accurately and conveniently than the original equations of motion. 

Using this approach Baumgarte et al. (1997d) studied the stability of binary neutron stars 
and found no evidence for the instability reported by Wilson and Mathews, even though they 
were using a very similar approximation scheme. One possible reason for the difference is that the 
models of Baumgarte et al. are corotating, that is, are locked in synchronous orbit. The models 
of Wilson and Mathews, by contrast, have very little intrinsic spin. It is possible in principle that 
the higher spin rate of the corotating models suppresses the instability ( [Mathews, Marronetti, 



fc Wilson 1997] ). As we shall see, the irrotational approximation allows the construction of 



models almost as easily as the corotating assumption does. However, the irrotational models are 
essentially nonspinning and should allow one to settle the collapse question definitively. 

Realistic neutron stars are in fact unlikely to be corotating — the viscosity is probably too low 



to synchronize the spin and orbital angular velocities ( Bildsten fe Cutler 1992 ; Kochanek 1992 



Lai 1994). Instead, the stars are likely to have very low intrinsic spins compared with the orbital 
frequency, and to good approximation we can assume they are irrotational (zero vorticity as seen 
from the inertial frame). Thus besides being of use in determining whether the Wilson-Mathews 
instability is real or not, the irrotational approximation will produce models that are close to those 
expected to actually occur. 

Ultimately, even if neutron star binaries do not collapse prematurely to black holes, they 
will reach an instability corresponding to the innermost stable circular orbit. The subsequent 
plunge and merger cannot be treated by the quasiequilibrium approximation. However, models 
constructed with the approximations in this paper should provide excellent initial data for the fully 
dynamical codes required to treat the final plunge. In fact, because of the extreme computational 
requirements of dynamical codes, it is unlikely that they will be able to follow the evolution from 
large separation at all. Realistic initial data such as provided by the methods here will be crucial. 



2. Irrotational Binaries in Newtonian Gravity 



Exact irrotational binaries in Newtonian gravity were first constructed by Bonazzola et al. 
(1992). They used the following method: Start with the equation of motion for a fluid in a 
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uniformly rotating frame of reference: 

<9v 1 

— + (v ■ V)v + 2ft x v + ft x (ft x r) = VP - V^. (1) 

ov po 

Here ft is the constant orbital angular velocity, P is the pressure, po the density, and <j> the 
Newtonian potential. The origin of the coordinates is at the center of mass of the binary system. 
The velocity u in the inertial frame is related to the velocity in the rotating frame by 

u = Hxr + v. (2) 

For later comparison with the relativistic case, note that the time derivatives in the rotating and 
inertial frames are related by 

— = - + £ w , (3) 

where w = ft x r is the orbital velocity and C denotes the Lie derivative. 

Now assume stationarity in the rotating frame: d/dt' — > 0. Then equation (|l|) can be 
rewritten as 

V(±v 2 + h N + 4> + 4> c ) + [(V x v) + 2ft] x v = 0. (4) 

Here 

r dP 

h N = / — = u + P po 5 
J Po 

is the Newtonian enthalpy per unit mass, u is the internal energy per unit mass, and 

1 ft x r) 2 (6) 



2 

is the centrifugal potential. Equation (|4|) can also be written in terms of the inertial velocity: 

V(i« 2 - (ft x r) • u + h N + (j)) + (V x u) x (u - ft x r) = 0. (7) 
The continuity equation, with d/dt' — > 0, becomes 

V • v = —v • , (8) 

Po 

or 

V-u = -(u-ft xr)--^. (9) 

Po 

The corotating case corresponds to v = 0. Equation (||) is trivially satisfied, while equation 
d|]) is just the Bernoulli equation: 

Jin + 4> + 4>c = const. (10) 

The irrotational case corresponds to no vorticity in the inertial frame, V x u = 0. This implies 
that u is given by a velocity potential, 

u = VV>. (11) 
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Equations ([?]) and (|9|) become 

\{V^f - (SI x r) • Vip + h N + (ft = const, (12) 

V 2 ^ = -(V^-fixr).^. (13) 

Po 

Equation (13) must be solved subject to the boundary condition 

(V^-ft x r) • VpoUrf = 0, (14) 

where the surface is defined by po = 0. 

Note that in terms of the velocity in the rotating frame, v, its divergence is given by equation 
(^) while its curl is 

V x v = -20. (15) 



The Bernoulli integral is 

2 



\v 2 + /in + (ft + (f) c = const. (16) 



The curl equation (15) is solved by 

v = -fixr + V^, (17) 

and we then recover equations (p^) and (13) from equations ( |Tfj| ) and (||). 

Equations ( |l2|) and (^) were solved by Bonazzola et al. (1992) for an illustrative case, but 
they did not give any sequences of models. The first sequences were constructed by Uryu & 
Eriguchi (1997), who solved the same equations for equal mass binaries with polytropic equations 
of state. 



3. Irrotational Flow in General Relativity 

Just as in Newtonian fluid mechanics, the velocity of a relativistic perfect fluid can be 



expressed as the gradient of a potential if the vorticity is zero (see, e.g., Landau & Lifshitz 195S, 
Moncricf 198C| ). Define the relativistic enthalpy by 



or 



h = P -±^, (18) 
Po 

where p is the total energy density and po = msn is the rest-mass density. (We use units with 
c = G = 1.) Here n is the baryon number density and mg the mean baryon rest mass. Then the 
relativistic vorticity tensor is defined as 

= P a ^u[^f3(hu a ) - V a {hu p )]. (19) 
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Here is the fluid 4-velocity and P^ u = 8^ v + u^u u is the projection tensor. There is another 
definition of relativistic vorticity, without the h in equation (|l^). Both definitions reduce to the 
correct Newtonian limit, where h — > 1, but the definition ( |l9| ) is the one that leads to a natural 
definition of potential flow. 

For a perfect fluid, Euler's equation can be written 

u Q V a (/iu^)+V M /i = 0. (20) 
Equations (|l^) and (^0|) yield a simple expression for the vorticity, 

= Vvihup) - VpQiUv). (21) 
Thus if the vorticity is zero, then the quantity hu u can be expressed as the gradient of a potential: 

hu^ = V^. (22) 

The equation of continuity for the baryon density n is 

V a (nu a ) = 0, (23) 

or, 

V Q [(n//i)V Q V] = 0. (24) 
The equation of state relates n to h, and h is found from the normalization condition u a u a = — 1, 



which by equation (22) gives 

/ l =[-(V Q ^)(V»] 1 / 2 . (25) 



Even for flow in a fixed background geometry, equation ( |24| ) with equation ( Pq) is in general 
a nonlinear equation in ijj and its derivatives. We rewrite it as 

V a V a ip = -(V>)V Q \n{n/h). (26) 



Then one way of solving equation (26) is by iteration, with either n/h or the whole right-hand 
side determined from the previous iteration. In general, this procedure will have to be iterated 
with the solution of the Einstein field equations for the metric. 



4. 3+1 Decomposition 

We make the usual ADM 3+1 decomposition of the spacetime metric: 
ds 2 = -a 2 dt 2 + j ij (dx i + f3 i dt){dx j + (3 j dt), 



(27) 
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where a is the lapse, f3 l is the shift, and 7™ is the 3-metric on the i = constant time slices. If we 
let n denote the unit normal to these time slices, and k the tangent to the t-coordinate lines (i.e., 
k = d/dt in the ADM coordinate system), the the decomposition is equivalent to writing 

k = an + p, n-0=O, (28) 

and expressing the 3-metric in a general coordinate system as 

Ivv = g^iv + n^riy. (29) 

We will call any object that is orthogonal to n spatial. 

The spatial (or 3-dimensional) covariant derivative of any tensor is defined by taking the 
4-dimensional covariant derivative and then projecting each free index with a so that the 
result is completely spatial. Since D^ a p = 0, this covariant derivative is compatible with y a p, 
and since a compatible covariant derivative is unique, this definition is equivalent to defining 
as the covariant derivative with respect to the 3-metric 7 Q/ g. 

For later reference, we list some standard relations that hold in the 3+1 decomposition. The 
derivative of the normal vector is 

V^n^ = -Kp V - n^Dy In a. (30) 
Here is the extrinsic curvature tensor, which is spatial. From equation (|3^) we see that 

nPV ll n v = D v ]jia. (31) 

The definition of gives 

V M / = IV - n^V u f, (32) 
where / is any scalar. Similarly, for any spatial vector X, we have 

V M X, = D^X V - n^n x V x X u - K^ x X x n u . (33) 



5. Quasiequilibrium 

As discussed in the Introduction, it is a good approximation to treat neutron star binaries 
as essentially equilibrium states, and to neglect the gravitational radiation on orbital time scales. 
The quasiequilibrium assumption is implemented mathematically by requiring that the spacetime 
admit a Killing vector to express the symmetry: a rotation A(f> about the rotation axis is equivalent 
to a displacement QAt, where £1 is the uniform orbital angular velocity. We write the Killing 
vector as 

f = — + n— = k + nm. (34) 
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Here the generator of rotations about the rotation axis of the binary is denoted by d/dcp, or fn in 
an arbitrary coordinate system. By equation (^), we can also write 

1 = an + B, (35) 

where 

B = f3 + nrh (36) 

is the rotating shift vector. It is equivalent to adding fi x r to the original shift. (Note that here B 
is the negative of B defined by ponazzola, Gourgoulhon, fc Marck 1997] .) While I is not a Killing 
vector in the exact metric, in the quasiequilibrium approximation Einstein's equations can be 
consistently approximated to neglect gravitational radiation and admit I as a Killing vector, e.g., 
by the Wilson-Mathews scheme flWilson fc Mathews 1989| ; [Wilson 1990) ; |Wilson fc Mathews 1995b - 



We require that I be a symmetry generator for the matter fields as well as for the metric. 
Mathematically, this means that the Lie derivative along I of any matter field must vanish (cf. eq. 
^). Note that since ip is a potential, we must be careful that the symmetry condition is expressed 
correctly for it. It is not true that 

If + n % = (wrong!) (37) 

Rather, 

= Cj{hu) = Cj{dip) = d{£pl>). (38) 

Here we have used equation (^) in the notation of differential forms, and used the fact that the 
Lie derivative of a form commutes with the exterior derivative. Thus as long as Cfp is a constant 
the symmetry will be respected. We write 

PVrf = ^ = ^ + ^ = -C, (39) 

where the minus sign is chosen for later convenience: C is a positive constant. 

Another way of seeing why equation ([39]) holds is to project the Euler equation (^) along I. 
Using Killing's equation 

+ V„Z M = 0, (40) 
and the fact that h satisfies the symmetry, we can derive the Bernoulli integral 

hu^ = constant. (41) 

By equation (p2[), this is equivalent to equation (|39|). 



One way of implementing equation (^) is to solve equation ( p6| ) by separating out the 
t-dependence with a solution of the form 



ip = -Ct + f(r,6,<f>-SH). 



(42) 
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We will proceed instead to derive a coordinate-independent equation that embodies the symmetry. 

We will need expressions for the normal derivatives of scalar quantities that satisfy the 
symmetry. For ip, equations (|35| ) and ( |39] ) give 

a 

= --(C + B^D^). (43) 

a 

For a scalar / that satisfies C^f = we have a similar equation with C = 0: 

= ——B^D^f. (44) 



6. 3+1 Decomposition of the Potential Equation 

We can now derive the 3+1 decomposition of the potential equation (|2^). First, equations 
and (|4|) give 

V ft il> = D lt i/> + n ft (C + B v D v if))/a. (45) 

Therefore 

V"V^ = VD^ + (V%)(C + B v B v f)ja + n"V M (<? + B v B v ^)ja. (46) 
Using equation (33) with X p — > B u ^, we can write the first term on the right-hand side of (|46| ) as 



rrDyft - n x V x {n»B^) + (n x V x n^)B^ 
D^D^ip - + (D^D* In a. 



(47) 

In the second term on the right-hand side of ( |4"6| ) we replace V^n^ by —if ( = —K^^). For the 

(48) 



third term, we use equation (44): 

n"V M -(C + B V B V ^) = --B^B^-iC + B u D 1/ ip). 
a a a 

Substituting equations (|47|) and ( [48| ) in equation (p6|), we find 

, , /»/•/> . • • ( ;> ,..;»/' i,. ,.. ^ 



V^V.V = B^Baip + {B,ip)B^\na (C + B u B u ip) 

Q 

5|) and a similar equation with ^ 



1 B^ 

or or 



(49) 



For the right-hand side of equation (26) we use equation 
replaced by hx[n/K) and C = 0. So the final form of equation p^ ) is 

D i + (D i i;)D i In a + (C + B^Djip) -^D { In a 

\ or n 



^Di(BW^) = - 



DV - (C + BW^) 



I?' 



L>, In 



(50) 
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Here we have used latin indices since the equation is completely spatial, with the metric being the 
3-metric 7^. For numerical work, it may be advantageous to rewrite equation (|5(]) in an equivalent 
form. For example, define 

A = C + BWjiP = a[h 2 + (AV0 2 ] 1/2 , (51) 
where the second equality follows from equation (|57]) below. Then equation (|(]) becomes 

D i D i ij} - ~ = - <D^ - A Bi ) A In (™) . (52) 



An equivalent form is obtained by the substitution 

aK = DiB\ (53) 
which follows from Killing's equation (^). Then 

DW^ - D~ = - (d^ - A B .) Di In (™\ . (54) 

Since a consistent approximation scheme for Einstein's equations with the assumptions made here, 
such as the Wilson-Mathews scheme, typically has K = 0, the form ( |52"|) may be better. 



At the surface of either star, n — > 0, h —> 1. Thus equation (|52) or ( |54j ) must be solved 
together with the boundary condition 



= 0. (55) 

surf 



a- 

The Bernoulli integral for the matter distribution follows from substituting equation ( [45] ) in 
equation (p5|): 

/> 2 = -(D^)D^ + A(c + BWjiP) 2 , (56) 
or 

or equivalently 

C = -BWrf + a[h 2 + (AV>) 2 ] 1/2 - (57) 
The positive sign for the square root in equation (|57]) is determined by the Newtonian limit (cf. 
§0)- 



7. Newtonian Limit 

It is easy to see that the above equations reduce to the expected form in the Newtonian limit. 
The Newtonian limit of the various quantities is 

Dit/j -» diij (58) 

B t -> (nxr)i (59) 

a -> 1 + (60) 

ft -» I + /17V (61) 

if 0. (62) 



Thus equation ( |57| ) becomes 




(63) 



In leading order, C — > 1. If we define C = 1 + C, then in next order 




(64) 




(65) 



8. Conclusions 



Equation (|5^), together with equation (|^) or ([57]), is the principal result of this paper. It 



is a relatively simple Poisson-like equation for a scalar field that determines the fluid velocity 
for irrotational flow in the quasiequilibrium approximation. Bonazzola, Gourgoulhon, & Marck 
(1997) have considered the same problem as treated here, but their formalism is considerably more 
complicated. Their results are compared with ours in Appendix A. We do note, however, that this 
paper was inspired by reading theirs. 

The equations derived here, together with a compatible approximation to Einstein's equations 
such as the Wilson-Mathews scheme, allow one to construct sequences of neutron star binaries that 
should be very good approximations to actual binaries during the inspiral phase. The formalism 
is not much more complicated than that used by Baumgarte et al. (1997a,b,c,d) to construct 
corotating sequences, since equations very similar to equation (|(]) are already being solved in that 
case. Work is underway to construct such more realistic sequences. 

After this paper was submitted for publication, I became aware of an independent paper 
by Shibata (1998) treating the same problem. His equations (2.18) and (2.22) are completely 
equivalent to our equations ([57]) and (|54|) when the appropriate substitutions are made. 



I thank Larry Kidder for helpful discussions. This work was supported in part by NSF Grant 
No. PHY 94-08378 to Cornell University and by the Grand Challenge Grant No. NSF PHY 
93-18152/ASC 93-18152. 
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A. Comparison with the Formalism of Bonazzola, Gourgoulhon, & Marck 

Bonazzola, Gourgoulhon, & Marck (1997) have developed a different way of treating 
irrotational quasiequilibrium binaries. Instead of working directly with the 4-velocity of the fluid 
as we have, they consider the 3-velocity in the rotating frame. They decompose the 4-velocity as 

u = T(v + V), (Al) 

where v oc I is the 4-velocity of the rotating frame, V is the 3-velocity in the rotating frame, and 
r = (1 — V ■ V)^ 1 / 2 is the Lorentz factor. Then they derive the Bernoulli integral Q4l"|), which in 
their variables is 

In h + $ + In T = const, (A2) 
where $ = ln(a 2 — B ■ B) 1 / 2 , and also equations for the divergence and curl of the velocity field, 

V-V = -V f ln(nr), (A3) 

V A V = -2lj + (V A V + 2lj) ■ V _ y A y$. (A4) 

V ■ V 

The "cross product" operation A used for the curl is defined in the rest frame of the rotating 
observer: 

(aAb) ll = v a e a ^yb\ (A5) 

The "rotation vector" uj is defined by 

V A v = 2oj. (A6) 
In the Newtonian limit, where to — > Q,, these equations reduce to equations (0), (|l5|), and (|l6|). 



Bonazzola, Gourgoulhon, and Marck carry out the 3+1 decomposition of the above equations 
to get 3-dimensional equations for the velocity field. They then decompose the velocity using 
both scalar and vector potentials, and obtain a scalar and a vector Poisson-like equation. In the 
Newtonian limit, the vector potential has an analytic solution corresponding to the term — ft x r 
in equation (|l7l) , leaving just the scalar potential to be solved for. However, there does not seem 
to be an analytic solution for the vector potential in the relativistic case, and so one has to solve 
for it numerically as well. 

The procedure of Bonazzola, Gourgoulhon, & Marck (1997) can actually be simplified 
somewhat by working with a slightly different 3-dimensional velocity p, defined by 

p = e~*V. (A7) 



Then equations (|A3| ) and ( |A4j ) become 

V-p = -V J? ln(nre*), (A8) 
VAp = -2e~®u. (A9) 



However, after the 3+1 decomposition and introduction of scalar and vector potentials, one is still 
left with complicated equations. 
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